home *** CD-ROM | disk | FTP | other *** search
/ NetNews Offline 2 / NetNews Offline Volume 2.iso / news / comp / lang / c++-part1 / 335 < prev    next >
Encoding:
Text File  |  1996-08-05  |  6.8 KB  |  193 lines
  1. Path: news.ioffe.rssi.ru!usenet
  2. From: "Eugene I Levin" <Levin@lh.ioffe.rssi.ru>
  3. Newsgroups: comp.lang.c,comp.lang.c++
  4. Subject: Re: RANDOM NUMBER GENERATOR
  5. Date: Thu,  4 Jan 1996 00:12:10 +0300 (MSK)
  6. Organization: Ioffe Institute
  7. Distribution: world
  8. Message-ID: <ABg4lwmmK4@lh.ioffe.rssi.ru>
  9. NNTP-Posting-Host: ccpti.ioffe.rssi.ru
  10. Mime-Version: 1.0
  11. Content-Transfer-Encoding: 8bit
  12. X-Mailer: dMail [Demos Mail for DOS v1.23]
  13.  
  14. BCNJ68B@prodigy.com (Tom Kellerman) wrote:
  15.  
  16. > I want to thank you all for your responses. I am going to try and put
  17. > something together. I have tried to use just the RAND() and SRAND()
  18. > functions , but for some reason the result I get after 1M+ rolls the
  19. > numbers I am getting do not appear random. Maybe I am wrong, I don't
  20. > know.
  21.  
  22. Tom, you are RIGHT! The congruent generator which is employed in
  23. rand() is the bad one! And everybody know that. I followed this
  24. discussion with a kind of growing surprise: people cite Knuth's book as
  25. if there was nothing new about random number generation during 25
  26. years! BTW, the book is old but very good, and the summary of 500+
  27. pages Knuth wrote about congruent generators is: "If you need a
  28. generator for something serious, use shuffling algorithm" ("Algorithm
  29. MIX" from his book, if I recall the name correctly). I mean Knuth's
  30. MIX, not the strange algorithm "from MS Developer Network" somebody
  31. cited. Knuth's MIX is good, but slow.
  32.  
  33. Knuth wrote his best in 1970, but then in late seventies the old shift-
  34. register algorithm was re-invented, and I thought that this closed the
  35. discussion about which generator is better (at least for simulating
  36. purposes: cryptography people have their own considerations). At least
  37. I use nothing but this algorithm since 1980 and it proved to be
  38. trustful in a few very complex multi-dimensional problems.
  39.  
  40. Here is source of C function I use (I also have a FORTRAN version).
  41. Use it and forget rand(). Further details you can find in "T.G.Lewis &
  42. W.H.Payne. Generalized Feedback Shift Register Pseudorandom Number
  43. Algorithm. - Journal of ACM, 1973, vol 20, N3, pp. 456-468.". Sorry, I
  44. don't have more modern references at hand, but you can find them in
  45. any university library computer searching for "Shift Register Random".
  46.  
  47. Best regards, Eugene
  48.  
  49. *--> CUT HERE <----> CUT HERE <----> CUT HERE <----> CUT HERE <--*
  50.  
  51. extern void  sr250s (unsigned);
  52. extern float sr250f ();
  53. extern int   sr250l (unsigned);
  54. extern int   sr250i (unsigned short);
  55.  
  56. /*
  57.  
  58. DESCRIPTION:
  59. Shift-register Random Number Generator.
  60.  
  61. Copyleft (c) 1987, 1992. Eugene I Levin, A.F.Ioffe institute,
  62.                          St. Petersburg, Russia.
  63.                          E-mail: Levin@lh.ioffe.rssi.ru
  64.  
  65. This is a public domain program and may be freely distributed under the
  66. condition of remaining the original Copyleft information.
  67.  
  68. The author would appreciate to be informed about any improvements of
  69. these functions, or their implementations for any other computer system.
  70.  
  71. TESTS:
  72. These functions have been tested on a i860 computer with GNU C compiler
  73. and on 80486 computer with Borland C++. The main generator engine and
  74. it's integer entries are more or less machine-independent. The float
  75. entry, sr250f() should work on any computer were the type float occupies
  76. 4 bytes and follows IEEE standards.
  77.  
  78. REALIZATION NOTES:
  79. This is a realization of standard shift-register algorithm for random
  80. bit streams:
  81.    x[n] = x[n-q] XOR x[n-p]                                        (1)
  82. where p and q satisfy the condition that the polynomial x^p + x^q + 1
  83. must be primitive on the GF(2) field. There are a few known pairs (p,
  84. q), p = 250 and q = 103 has been used.
  85.  
  86. In order to generate float numbers independent bit streams (1) are used
  87. for the each bit of mantissa.
  88.  
  89. Consecutive values of (1) are stored in array "buffer" of the length p.
  90. The index of the n-th member of (1) in buffer is
  91.    i = n mod p                                                     (2)
  92. so that a new value x[n] always replaces x[n-p]. Substitution of Eq. (2)
  93. into Eq. (1) yields to the algorithm:
  94.    x[i] = x[i] XOR x[i - p]      when i >= p,                      (3a)
  95.    x[i] = x[i] XOR x[i - p + q]  when i < p.                       (3b)
  96.  
  97.                                                                    */
  98. #include <stdlib.h>
  99.  
  100. #define P  250
  101. #define Q  103
  102.  
  103. #define Q1             (P-Q)
  104. #define MANTISSA_MASK  0x003fffffl
  105. #define INT_MASK       0x0000ffffl
  106. #define FLOAT_2        0x40000000l  /* (float) 2  */
  107.  
  108. static union {
  109.              float f;
  110.              long  l;
  111.              }  mixer;
  112.  
  113. static long buffer[P];   /* the buffer described above */
  114. static int  index;       /* pointer to the current cell in the buffer */
  115. static int  i;           /* working variable */
  116.  
  117. /*
  118.    The following macro is the engine of the generator. It tests whether
  119.    there is an unused value in the buffer and, if there is not, replaces
  120.    the whole buffer with fresh values according to formulae (3).
  121.  
  122. */
  123.  
  124. #define REFRESH                                      \
  125.   if (index == P)                                    \
  126.      {                                               \
  127.      index = 0;                                      \
  128.      for (i=0; i<Q; i++) buffer[i] ^= buffer[i+Q1];  \
  129.      for (i=Q; i<P; i++) buffer[i] ^= buffer[i-Q];   \
  130.      }
  131.  
  132. /*====================================================================*/
  133.  
  134. void sr250s(unsigned seed)
  135.  
  136. /* This function fills the buffer with initial values taken from the
  137.    system generator (members -p through 0 of the sequence (1)).
  138.    sr250u *MUST* be called before any call to other sr250 functions.
  139.    Parameter "seed" determines the sequence of random numbers. */
  140.  
  141. {
  142. #define longrand ((long) rand())
  143.  
  144.   srand(seed);
  145.   for (i=0; i<P; i++)
  146.      buffer[i] = (longrand | (longrand<<15)) & MANTISSA_MASK;
  147.   index = P;  /* Marks the buffer as already used */
  148.  
  149. }
  150.  
  151. /*====================================================================*/
  152.  
  153. float sr250f()
  154. /* Generates a random float 0 <= x < 1. */
  155. {
  156.    REFRESH
  157.    mixer.l = buffer[index++] | FLOAT_2;
  158.    return (mixer.f - 2.0);
  159.    /* Floats on computers which follow IEEE float-point standards
  160.       (Intel, Sun) are stored without first bit in their mantissas
  161.       (which is supposed to equal 1), thus mixer.f belongs to the
  162.       interval [2, 3], not to [2, 4] as it would be on a "classic"
  163.       computer.    */
  164. }
  165.  
  166. /*====================================================================*/
  167.  
  168. int sr250i(unsigned short m)
  169.  
  170. /* Generates an integer 0 <= i < m. */
  171.  
  172. {
  173.    REFRESH
  174.    return (int) (((buffer[index++] & INT_MASK)*m)>>16);
  175. }
  176.  
  177. /*====================================================================*/
  178.  
  179. int sr250l(unsigned l1)
  180. /*
  181.    Generates an integer 0 <= i <= l1, where l1 *MUST* equal 2^n - 1.
  182.    This entry is somewhat faster than sr250i.
  183. */
  184. {
  185.    REFRESH
  186.    return (int) (buffer[index++] &((long) l1));
  187. }
  188.  
  189. --
  190. ===============================================================
  191. Eugene I Levin, St.Petersburg, Russia (Levin@lh.ioffe.rssi.ru).
  192.  
  193.